TUTORIAL 6 : Free Energy (FE) and replica-exchange methods

Author:Andrey V. Brukhno (andrey.brukhno{at}stfc.ac.uk)

General concept

Most often it is impossible to calculate the absolute free energy value (we will see why later on). However, relative free energy can always be expressed via the probability of the associated state:

\[{\rm Free\ energy\ of\ state}\ = \ -kT \ln{\rm (probability\ of\ state)} \ \ \{+\ unknown\ const\}\]

then for the FE difference we get

\[{\rm FED}_{12} = -kT \ln\frac{\rm probability\ of\ state\ 2}{\rm probability\ of\ state\ 1}\]

This is the most useful statistical-mechanical approach to FE evaluation, and all advanced FE methods are based on the above formula (althogh there exist also thermodynamic approaches, usually more demanding and time-consuming).

Example: Potential of mean force as FED

Potential of mean force (PMF or POMF) is defined as the work done upon reversibly bringing two particles from the state of zero-interaction into the state of full-interaction at a given distance, W(R).

PMF can be directly obtained from RDF:

\[\beta W(R) = -\ln g(R) + \ln g(\infty)\]

where the zero-interaction state corresponds to infinite separation \((R\to\infty)\).

Alternatively, one can calculte

\[\beta \Delta W(R_1 \to R_2) = -\ln g(R_2) + \ln g(R_1) = -\ln \frac{g(R_2)}{g(R_1)}\]
  • the work to bring two particles \((i,j)\) from \(|\Delta\mathbf r_{ij}|=R_1\) to \(R_2\).

Clearly, if \(g(R_1)=1\), we get back the total work \(W(R=\infty\to R_2)\).


Why PMF is a free energy?

First, note that

\[\frac{g(R_2)}{g(R_1)} = \frac{ \langle \delta(|\Delta\mathbf r_{ij}|-R_2) \rangle } { \langle \delta(|\Delta\mathbf r_{ij}|-R_1) \rangle } = \frac{ \int \delta(|\Delta\mathbf r_{ij}|-R_2)\exp[-\beta U(\mathbf r^{N})] \ d{\mathbf r}^{N} } { \int \delta(|\Delta\mathbf r_{ij}|-R_1)\exp[-\beta U(\mathbf r^{N})] \ d{\mathbf r}^{N} }\]

We can rewrite the canonical configuration integral (configurational part of the partition function)

\[Q(N,V,T) = \int dR \int \delta(|\Delta\mathbf r|-R)\exp[-\beta U(\mathbf r^{N})] \ d{\mathbf r}^{N} = \int Q_{NVT}(R) dR\]

as an integral (i.e. sum) over the substates corresponding to \(|\Delta\mathbf r_{ij}|=R\). This implies that we can consider each substate as an individual thermodynamic state, whereas \(R\) serves as an additional external constraint, similar to N,V,T or p. Thereby,

\[\beta \Delta W(R_1 \to R_2) = -\ln \frac{g(R_2)}{g(R_1)} = \frac{ \langle \delta(|\Delta\mathbf r_{ij}|-R_2) \rangle } { \langle \delta(|\Delta\mathbf r_{ij}|-R_1) \rangle } = -\ln \frac{Q_{NVT}(R_2)}{Q_{NVT}(R_1)} = -\ln \frac{p(R_2)}{p(R_1)}\]


\[p(R) = \frac{Q_{NVT}(R)}{Q(N,V,T)}\]

Therefore, we can evaluate the FED between any two substates (different particle separations) by performing a Metropolis sampling along R and collecting a histogram of hits (or visits) on a discrete grid with sufficiently small bin size, \(\Delta R\), which is of course equivalent to calculating \(g(R)\)!

NOTE: One does not need to know the total partition function (nor the total configuration integral) as it only serves as a normalization constant which cancels out!

  • Question: what is the “unknown constant” in the first equation above?
  • Similar formulae hold true generally, for other constrained substates (see the Appendix below).

FED methods in DL_MONTE

In general, Metropolis sampling is prone to getting trapped (stuck) in local FE minima and, hence, suffering from undersampling (even avoiding!) regions around energetic or entropic barriers. This makes estimating FED via probablity ratio troublesome…

The solution to this problem is to bias the sampling away from the minima and towards the barriers.

  • Ideally the bias potential is to perfectly compensate for the underlying free energy differences!
  • This implies that the resulting probablity distribution in a perfectly biased simulation is uniform, i.e. flat!

Parallel tempering with replica-exchange (RE)

Two approaches to parallelization in DL_MONTE


In a parallel multicanonical simualtion each workgroup treats the system at a different tempreature and generates it own MC trajectory. It is then possible to arrange periodic exchnage of configurations between the workgroups, which is equivalent to concurrent jumps in T with the acceptance probability:

\[acc[ \{ \mathbf{r}^N; T_1 \} \leftrightarrow \{ \tilde{\mathbf r}^N; T_2 \} ] = \min(1, \exp \{- (\beta_1 U_{2} - \beta_2 U_2) - (\beta_2 U_1 - \beta_1 U_1) ] \} ) = \min(1, \exp \{- \Delta \beta_{12} \Delta U_{12} \} )\]

If a FED evaluation is carried out in the same simulation, say FED(R), the acceptance rule includes an extra term:

\[acc[ \{ \mathbf{r}^N; T_1 \} \leftrightarrow \{ \tilde{\mathbf r}^N; T_2 \} ] = \min(1, \exp \{- \Delta \beta_{12} \Delta U_{12} - \Delta \beta_{12} \Delta U_{12}^{(bias)} \} )\]


Now try the exercises to learn how to perform FED calculations with DL_MONTE:

TUTORIAL 6 : Exercise 1 - Reproducing interaction potentials U(r) with EE & WL/RE methods

TUTORIAL 6 : Exercise 2 - Calculating PMF for two charged colloids in ionic cloud

TUTORIAL 6 : Exercise 3 - Umbrella sampling in subranges (windows) with the use of WHAM method

Appendix: FED via a Generalised (Expanded) Ensemble

The biasing is best illustrated with a generalised, so-called Expanded Ensemble (EE) approach.

The EE is defined by a generalised partition sum:

\[Z_{EE} = \int w({\mathbf q}) d{\mathbf q} \int_{ens} \exp[-\beta H({\mathbf r}^N; {\mathbf q})] d{\mathbf r}^N\]

or in a discretised form:

\[Z_{EE} = \sum_m w({\mathbf q_m}) \int_{ens} \exp[-\beta H({\mathbf r}^N; {\mathbf q_m})] d{\mathbf r}^N\]

where \(\{{\mathbf q}\}\) stand for one or more generalised coordinates (or constraints), which can also include normal thermodynamic variables (e.g. N,V,T, p etc), and \(w({\mathbf q})\) is a biasing weight assigned to a particular set of \({\mathbf q}\).

Clearly, if Metropolis sampling is carried out in such a generalised (so-called expanded) ensemble, including MC moves in \({\mathbf q}\), the probablity of an EE substate

\[p({\mathbf q}) = \frac{w({\mathbf q}) d{\mathbf q} \int_{ens} \exp[-\beta H({\mathbf r}^N; {\mathbf q})] d{\mathbf r}^N}{Z_{EE}} = \frac{w({\mathbf q}_m) \int_{ens} \exp[-\beta H({\mathbf r}^N; {\mathbf q}_m)] d{\mathbf r}^N}{Z_{EE}}\]

and the effective (biased) free energy difference

\[\beta \Delta F_{nm}^{(EE)} = -\ln \frac{p({\mathbf q}_m)}{p({\mathbf q}_n)} = -\ln \frac{w({\mathbf q}_m) Q_{ens}({\mathbf q}_m)}{w({\mathbf q}_n) Q_{ens}({\mathbf q}_n)} = -\ln \frac{w({\mathbf q}_m)}{w({\mathbf q}_n)} + \beta \Delta F_{nm}^{(true)}\]

Finally, if the weights \(\{w({\mathbf q}_m)\}\) are chosen so that \(p({\mathbf q}_m)/p({\mathbf q}_n) = 1\), we get:

\[\beta\Delta F_{nm}^{(true)} - \ln \frac{w({\mathbf q}_m)}{w({\mathbf q}_n)} = 0\]
  • Ideally the EE weights are to perfectly compensate for the unbiased probablity distribution!

In practice, though, one has to arrange an iteration to self-consistently update the weights until all the substates are populated reasonably well, after which the FED is obtained as

\[\beta\Delta F_{nm}^{(true)} = \ln \frac{w({\mathbf q}_m)}{w({\mathbf q}_n)} - \ln \frac{p({\mathbf q}_m)}{p({\mathbf q}_n)} = \ln \frac{w({\mathbf q}_m) p({\mathbf q}_n)}{w({\mathbf q}_n)p({\mathbf q}_m)}\]

NOTE: It is convient (and conventional) to define the EE weights via a biasing potential:

\[w({\mathbf q}) = \exp[-\beta U^{(bias)}({\mathbf q})]\]

so that

\[\beta\Delta F_{nm} = -\Delta U_{nm}^{(bias)}({\mathbf q}) - \ln \frac{p({\mathbf q}_m)}{p({\mathbf q}_n)}\]
  • Ideally the bias potential is to perfectly compensate for the underlying free energy differences!